Meron-Cluster Solution of Fermion Sign Problems 
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We present a general strategy to solve the notorious fermion sign problem using cluster algorithms. 
The method applies to various systems in the Hubbard model family as well as to relativistic 
fermions. Here it is illustrated for non-relativistic lattice fermions. A configuration of fermion 
world-lines is decomposed into clusters that contribute independently to the fermion permutation 
sign. A cluster whose fiip changes the sign is referred to as a meron. Configurations containing 
meron-clusters contribute to the path integral, while all other configurations contribute 1. The 
cluster representation describes the partition function as a gas of clusters in the zero-meron sector. 



CD 






I 

o 
o 



> 

00 

O 



X3 
O 

o 



X 



02.70.Lq, 71.10.Fd, 05.50+q, 12.38.Gc 

The numerical simulation of fermions is a notorious 
problem that hinders progress in understanding high- 
temperature superconductivity [^], QCD at non-zero 
chemical potential Ig] and many other important prob- 
lems in physics. One of the main problems originates 
from the minus-signs associated with Fermi statistics 
which prevent us from interpreting the Boltzmann fac- 
tor in a fermionic path integral as a positive probability. 
When the sign of the Boltzmann factor is incorporated 
in measured observables, the fluctuations in the sign give 
rise to dramatic cancellations. Especially for large sys- 
tems at low temperatures this leads to relative statistical 
errors that are exponentially large in both the volume 
and the inverse temperature. This makes it impossible in 
practice to study such systems with standard numerical 
methods. Here, for the first time, we completely elim- 
inate a severe sign problem in the simulation of a non- 
relativistic system of interacting lattice fermions using a 
cluster algorithm. The solution of the problem proceeds 
in two steps. The idea of the first step is to use cluster 
algorithm techniques to reduce the problem of canceling 
many contributions ±1 to the problem of averaging over 
non-negative contributions and 1. This step solves one 
half of the sign problem as discussed below. In large 
volumes and at small temperatures one still generates 
vanishing contributions to the average sign most of the 
time and very rarely one encounters a contribution 1. In 
order to solve the other half of the problem a second step 
is necessary which guarantees that contributions and 1 
are generated with similar probabilities. The idea behind 
the second step is to include a Metropolis decision in the 
process of cluster decomposition. The two basic ideas 
behind our algorithm are general and apply to a variety 
of systems. In this paper, we illustrate our method for 
a simple model which serves as a testing ground for the 
new ideas. 

Let us consider a fermionic path integral Zf = 
^„ Sign[ri] exp(— S'[n]) over configurations n with a 
Boltzmann weight of Sign[?i] = ±1 and magnitude 
exp(— S'[n]). Here S[n] is the action of a corre- 



sponding bosonic model with partition function Zi, = 
X]„exp(— S'[n]). A fermionic observable 0[n] is obtained 
in a simulation of the bosonic ensemble as 
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^^0[n]Sign[n]exp(-5M) = ^^^. (1) 



The average sign in the simulated bosonic ensemble is 



(Sign) = |t = cxp{-f3VAf). 
Zb 



(2) 



The last equality (valid for large (3V) points to the heart 
of the sign problem. The expectation value of the sign 
is exponentially small in both the volume V and the in- 
verse temperature /3 because the difference between the 
free energy densities A/ = ff — ft of the fermionic and 
bosonic systems is necessarily positive. Even in an ideal 
simulation of the bosonic ensemble which generates N 
completely uncorrelated configurations, the relative sta- 
tistical error of the sign (again for large (3V) is 



ASign _ yiS^g^ > - (Sign)2 ^ cxp{f3VAf) 
(Sign) 



iV(Sign) 



(3) 



Here we have used Sign =1. To determine the aver- 
age sign with sufficient accuracy one needs to generate 
on the order of A^ = exp{2f3VAf) configurations. For 
large volumes and small temperatures this is impossible 
in practice. 

It is possible to solve one half of the problem if one 
can match all contributions —1 with 1 to give 0, such 
that only a few unmatched contributions 1 remain. Then 
effectively Sign = 0, 1 and hence Sign — Sign. This 
reduces the relative error to 



ASign _ V(Sign) - (Sign)^ _ exp(/3yA//2) 



(Sign) 



/]V'(Sign) 
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One gains an exponential factor in statistics, but one 
still needs to generate A^' — \/N = exp(/?yA/) indepen- 
dent configurations in order to accurately determine the 



average sign g. This is because one generates exponen- 
tially many vanishing contributions before one encoun- 
ters a contribution 1. 

In several cases cluster algorithms provide an explicit 
matching of contributions —1 and 1 using an improved 
estimator. Cluster algorithms are a very efficient tool to 
simulate quantum spin systems [0-p|- In particular, the 
method can be implemented directly in the Euclidean 
time continuum [Q. The basic idea behind these algo- 
rithms is to decompose a configuration into Nc clusters 
of spins which can be flipped independently. Averaging 
analytically over the 2^^ configurations generated by the 
cluster flips, one can construct improved estimators for 
various physical quantities. As we will show, using an im- 
proved estimator for the fermion sign, cluster algorithms 
can solve the sign problem if the clusters contribute inde- 
pendently to the sign and a reference cluster orientation 
with a positive weight always exists. This means that the 
flip of any given cluster either changes the sign or not, 
independent of the orientation of all the other clusters. A 
cluster algorithm for lattice fermions was first presented 
in PI with the hope of finding such an improved esti- 
mator. Unfortunately, in that algorithm the clusters do 
not affect the sign independent of one another. Still, 
cluster algorithms have been used for fermion models P| . 
For systems with no severe sign problem these algorithms 
work much better than standard numerical methods, but 
they do not solve the fermion sign problem. 

A solution to a sign problem using cluster algorithms 
was first found in a bosonic model with a complex ac- 
tion — the 2-d 0(3) model at vacuum angle 9 — n [nOJ. 
The cluster independence was achieved by constructing 
a non-standard action. In that model clusters whose 
flip changes the sign are half-instantons which are usu- 
ally referred to as merons. In this paper we extend the 
meron-concept to fermionic models by demanding clus- 
ter independence. For non-relativistic spinless fermions 
hopping on a d-dimensional cubic lattice of size V = L''' 
with periodic boundary conditions, this leads us to the 
Hamiltonian lllll 
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with [/ > i > 0. Here i is a unit vector in the i-direction, 
c^ and Cx are fermion creation and annihilation opera- 
tors obeying the standard anticommutation relations and 
Ur^ — c^Cx is the occupation number of the lattice site 
X. Since f/ > 0, two fermions or two holes on neighbor- 
ing lattice sites repel each other, while a fermion and a 
hole attract one another. This is a simple example of a 
fermionic model for which the sign problem can be solved 
completely using a meron-cluster algorithm. 

Let us now discuss our model and algorithm in more 
detail. Following S we introduce a space-time lattice 



with 2dM time-slices and spacing e = (3/M in the Eu- 
clidean time direction and we insert complete sets of oc- 
cupation number n(x, t) = 0, 1 eigenstates at each time- 
slice to express the partition function as a path integral. 
The magnitude exp(— 5'[n]) of the Boltzmann factor is a 
product of four-site interactions associated with space- 
time plaquettc configurations [n(x^ t), n{x -I- z, t), n{x, t + 
l),n(a; + i,t + 1)]. The sign factor Sign[ri] — ±1 has a 
topological meaning. The occupied sites form fermion 
world-lines which are closed in Euclidean time. Parti- 
cles may be exchanged during their Euclidean time evo- 
lution and the fermion world-lines define a permutation 
of particles. According to the Pauli principle, Sign[n] is 
just the sign of that permutation. In the following we 
restrict ourselves to U — t. Then the bosonic system 
without the sign factor is the antiferromagnetic spin 1/2 
quantum Heisenberg model, and n{x^ t) — Q and 1 cor- 
respond to spin —1/2 and 1/2, respectively. The stag- 
gered occupation (the analog of the staggered magneti- 
zation) 0[n] = eX;^ j(-l)^i+^2+...-h:rd(^(a.^i) _ i)^ and 

the corresponding susceptibility x = (O^ Sign) //3T^ (Sign) 
are important observables. 

The algorithm decomposes a configuration into closed 
loops of lattice points which may be flipped indepen- 
dently. When a loop is flipped, the occupation numbers 
of all points on the loop are changed from to 1 and vice 
versa. Each lattice point participates in two space-time 
plaquette interactions [n{x,t), n(x-\-i,t), n{x,t+l), n(x+ 
i,t+ 1)]. On each interaction plaquette the lattice points 
are connected in pairs and a sequence of connected points 
defines a loop-cluster. For space-time plaquette configu- 
rations [0, 0, 0, 0] and [1, 1, 1, 1] the lattice points are con- 
nected with their time-like neighbors, for configurations 
[0,1,1,0] and [1,0,0,1] they are connected with their 
space-like neighbors and for configurations [0, 1, 0, 1] and 
[1,0, 1,0] they are connected with their time-like neigh- 
bors with probability p = 2/(1 + exp(eC/)) and with their 
space-like neighbors with probability I — p. After iden- 
tifying the clusters, they are flipped independently with 
probability 1/2. 

A remarkable property of the cluster rules is that 
Sign[n] = riiJi Sign[Ci], where Ci,i — l,...,Nc denotes 
the oriented clusters in a conflguration. By properly flip- 
ping the clusters, one can reach a reference configuration 
(the first configuration in figure 1) in which all even lat- 
tice sites are occupied and all odd sites are empty. In the 
reference orientation Sign[Ci] = 1. When the cluster is 
flipped, Sign[Ci] = 1 if N^, + Nii/2 is odd and —1 other- 
wise. Here iV„j is the temporal cluster winding number 
and Nh is the number of times the cluster hops to a neigh- 
boring lattice point. This relation follows directly from 
the fermionic anticommutation relations. Following [[lO| , 
we refer to clusters whose flip changes the sign as merons. 
The flip of a meron-cluster permutes the fermions and 
changes the topology of the fermion world-lines. Since 



flipping all clusters does not change the fermion sign, the 
number of meron-clusters is always even. Two fermion 
configurations together with a meron-cluster are shown 
in figure 1. 

Sign[n] = f 



X- 



{Y.c\Oc?5Nfi + 2\OcA\Oc.\5^,2) 
VP{5n.») 



(6) 




X 

Sign[n] = -f 



1 i 


\ 


< 




Ik 


1 


( 

k 







FIG. 1. Two configurations of fermion occupation numbers 
in (1 + 1) dimensions. The shaded plaquettes carry the inter- 
action. The dots represent occupied sites. In the second con- 
figuration two fermions interchange their positions. Flipping 
a meron-cluster (represented by the fat line) changes one con- 
figuration into the other and changes the fermion sign. The 
other clusters in the configurations are not shown. 

The improved estimator for (Sign) is the average over 
the 2 '^ configurations obtained from independently flip- 
ping the Nc clusters in all possible ways. For conflgu- 
rations that contain merons the average sign is zero be- 
cause flipping a single meron leads to a cancellation of 
signs ±1. Only the configurations without merons con- 
tribute to (Sign) and their contribution is always 1. This 
solves one half of the sign problem as discussed before. 

Let us now consider an improved estimator for 
(O^^Sign) which is needed to determine the susceptibil- 
ity X- The staggered occupation, 0[n] = J2c^c, is 
a sum of staggered occupations of the clusters, Oc — 

eE(x.t)ec(-l)'''^'''^'"^''''("(2;'*) - l)- When a cluster 
is fiipped, its staggered occupation changes sign. In a 
configuration without merons, where Sign[n] = 1 for all 
relative cluster flips, the average of 0[n]^Sign[n] over all 
2^'^ configurations is ^^ lOcp. For configurations with 
two merons the average is 210,^^110(72 1 where Oi and O2 
are the two meron-clusters. Configurations with more 
than two merons do not contribute to (O^Sign). Thus, 
the improved estimator for the susceptibility is given by 



where N is the number of meron-clusters in a configura- 
tion. Hence, to determine x 0116 must sample the zero- 
and two-meron sectors only. 

The probability to find a configuration without merons 
is exponentially small in the space-time volume since it 
is equal to (Sign). Thus, although we have increased 
the statistics tremendously with the improved estimators, 
without a second step one would still need an exponen- 
tially large statistics to accurately determine x- One goal 
of the second step is to eliminate all configurations with 
more than two merons. This enhances both the numer- 
ator and the denominator in eq.(|6|) by an exponentially 
large factor, but leaves their ratio unchanged. We start 
with an initial configuration with zero or two merons. 
For example, a completely occupied configuration has no 
merons. We then visit all plaquette interactions one af- 
ter the other and choose new cluster connections between 
the four sites according to the cluster rules. If the new 
connection increases the number of merons beyond two, 
it is not accepted and the old connection is kept for that 
plaquette. To decide if the meron number changes, one 
needs to examine the clusters affected by the new con- 
nection. Although this requires a computational effort 
proportional to the cluster size (and hence to the phys- 
ical correlation length) this is no problem, because one 
gains a factor that is exponentially large in the volume. 
The above procedure obeys detailed balance because con- 
figurations with more than two merons do not contribute 
to the observables we consider. Also, one can show that 
the algorithm is still ergodic. The simple reject step elim- 
inates almost all configurations with weight and is the 
essential step to solve the other half of the fermion sign 
problem. 

Since for large space-time volumes the two-meron sec- 
tor is much larger than the zero-meron sector, with- 
out further improvements one would still need statistics 
quadratic (but no longer exponential) in the space-time 
volume to accurately measure x- The remaining problem 
can be solved with a re- weighting technique similar to the 
one used in pO| . To enhance the zero-meron configura- 
tions in a controlled way, we introduce a trial probability 
Pt{N) for each iV- meron sector. We set pt{N) for N > 2 
to infinity and use it in a Metropolis accept-reject step 
for the newly proposed cluster connection on a specific 
plaquette interaction. A new connection that changes 
the meron number from N to A'"' is accepted with prob- 
ability p = nim[l,pt{N)/pt{N')]. In particular, config- 
urations with N' > 2 are never generated because then 
Pt{N') — 00 and p — 0. After visiting all plaquette in- 
teractions, each cluster is flipped with probabihty 1/2 
which completes one update sweep. After re-weighting, 
the zero- and two-meron configurations appear with sim- 
ilar probabilities. This completes the second step in our 
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TABLE I. Numerical results for {Sign) and x obtained with 
algorithm A^ and x obtained with algorithm A2. 



solution of the fermion sign problem. The re-weighting 
of the zero- and two-meron configurations is taken into 
account in the final expression for the susceptibility as 
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We have implemented the meron cluster algorithm in 
(24-1) dimensions and have tested it using exact diagonal- 
ization results on small lattices. Table 1 contains a com- 
parison of results obtained with two algorithms using the 
same number of sweeps in both cases. The first algorithm 
(Ai) has the improved estimators and solves one half of 
the sign problem. The second algorithm (A2) has both 
the improved estimators and the additional Metropolis 
step and also solves the other half of the problem. The 
algorithm A2 is clearly superior once the average sign 
becomes small. In particular, we have applied A2 to sys- 
tems of size V — 12^ at a low temperature (3U — 8. This 
is far beyond reach of standard fermion algorithms and 
even of the algorithm Ai. It should be noted that our 
model has a very severe sign problem which persists after 
integrating out the fermions even at half-filling. 

Cluster representations in general and the meron- 
concept in particular are more than mere algorithmic 
tools. In fact, we have shown that the fermionic par- 
tition function can be expressed as a classical statistical 
mechanics system of clusters. The cluster formulation is 
a novel type of bosonization which works in any dimen- 
sion. In this formulation the Pauli principle manifests 
itself by the vanishing Boltzmann weight of a configura- 
tion containing meron-clusters. If we ignore the fermion 
permutation sign, the theory describes a gas of merons 
and non-merons with a large configuration space. Includ- 
ing the sign factor forces even numbers of merons to be 
bound into non-merons. As a consequence, in agreement 
with the Pauli principle, the configuration space is very 
restricted. The merons allow us to simulate fermions with 
local bosonic variables. This is much more efficient than 
integrating out the fermions, which leads to non-local 
bosonic effective actions. 

While the details of our algorithm are specific to the 
fermion model we have considered, the two basic ideas be- 
hind it are general and apply to a variety of models. They 



lead to a complete solution of the fermion sign problem 
for models of relativistic staggered fermions [^ as well as 
for non-relativistic fermions with spin. In applications of 
the meron-cluster algorithm to systems in the Hubbard 
model family we have so far not found high-temperature 
superconductivity. Meron-cluster algorithms are also ap- 
plicable to quantum spin models in an arbitrary mag- 
netic field for which a similar type of sign problem arises. 
Similarly, one can solve the sign problem resulting from a 
complex action in the 2-d 0(3) model at non-zero chem- 
ical potential or at non-zero vacuum angle 0. The next 
challenge is to find applications of this method to QCD 
at non-zero baryon density. It seems likely that progress 
along the lines discussed here can be made in the quan- 
tum link D-theory formulation of the problem [|l2|,^ . 

U.-J. W. likes to thank the physics department of Duke 
University where part of this work was done for its hospi- 
tality and the A. P. Sloan foundation for its support. This 
work is supported in part by funds provided by the U.S. 
Department of Energy (D.O.E.) under cooperative re- 
search agreements DE-FC02-94ER40818 and DE-FG02- 
96ER40945. 



[1] S. R. White et al., Phys. Rev. 40 (1989) 506; 

E. Dagotto et al, Phys. Rev. B41 (1990) 811. 
[2] I. Barbour et al., Nucl. Phys. B (Proc. Suppl.) 60A (1998) 

220; 

M. Alford, A. Kapustin and F. Wilczek, Phys. Rev. D59 

(1999) 054502. 
[3] The fact that an improved estimator alone cannot solve 

the sign problem was pointed out to one of the authors 

by H. G. Evertz a long time ago. 
[4] U.-J. Wiese and H.-P Ying, Phys. Lett. A168 (1992) 143. 
[5] H. G. Evertz, G. Lana and M. Marcu, Phys. Rev. Lett. 

70 (1993) 875. 
[6] U.-J. Wiese and H.-P. Ying, Z. Phys. B93 (1994) 147. 
[7] B. B. Beard and U.-J. Wiese, Phys. Rev. Lett. 77 (1996) 

5130. 
[8] U.-J. Wiese, Phys. Lett. B311 (1993) 235. 
[9] N. Kawashima, J. E. Gubernatis and H. G. Evertz, Phys. 

Rev. B50 (1994) 136; 

B. Ammon et al., Phys. Rev. B58 (1998) 4304; 

M. Brunrier and A. Muramatsu, Phys. Rev. B58 (1998) 

RIOIOO. 
[10] W. Bietenholz, A. Pochinsky and U.-J. Wiese, Phys. Rev. 

Lett. 75 (1995) 4524. 
[11] S. Chandrasekha ran, J. Cox, K. Holland and U.-J. Wiese, 

[12 

[13 



hep-lat/9906021 



S. Chandrasekharan and U.-J. Wiese, Nucl. Phys. B492 
(1997) 455. 

R. Brower, S. Chandrasekharan and U.-J. Wiese, pep 
th/9704106l to appear in Phys. Rev. D. 



